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We propose a lattice model for Dirac fermions which allows us to break the degeneracy of the 
node structure. In the presence of a random gap we analyze the scaling behavior of the localization 
length as a function of the system width within a numerical transfer-matrix approach. Depending 
on the strength of the randomness, there are different scaling regimes for weak, intermediate and 
s ! . strong disorder. These regimes are separated by transitions that are characterized by one-parameter 

scaling. 
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I. INTRODUCTION 



Two-dimensional Dirac fermions play a crucial role in graphene and on the surface of topological insulators. A 
fascinating observation in graphene is the robust electronic transport in the vicinity of the two Dirac nodes, where 
two electronic bands meet each other with linear dispersion [l], Q. The latter is a consequence of the honeycomb 
lattice in graphene, which decomposes into two triangular lattices. 

In contrast to the experimentally observed robust transport properties it has been claimed from the theoretical 
side that transport is very sensitive whether inter-node scattering is present or not in the presence of disorder Q. In 
^ \ particular, there has been speculations that electronic states are delocalized in the absence of inter-node scattering but 
) • . localized in its presence. This has been explained by changing the symmetry class of the underlying Hamiltonian from 
symplectic to orthogonal [1, Q • These claims are based on weak- localization calculations [!, Q , which predict weak 
(anti-) localization (with) without inter-node scattering. Since weak localization calculations can only indicate the 
tendency towards localization, it would be interesting to evaluate this effect directly in terms of the scaling behavior 
of the localization length. For this purpose we shall study the localization length of a strip of finite width M under 
a change of M in this paper. Our method, originally introduced for transfer-matrix calculations of the Schrddinger 
Hamiltonian @,0|, will be applied subsequently to 2D lattice Dirac fermions with one or more nodes. For this purpose 
we introduce a model which has two bands and four Dirac nodes. We can open a gap at one node and gaps for the 
other three nods independently. This allows us to study the effect of intervally scattering by either keeping all four 
nodes or removing three of them and keeping only a single node. 
■ The aim of this work is to understand the scaling behavior of the localization length in two dimensions in the 
, metallic regime and near a metal-insulator transition due to a gap opening. The latter has been observed recently 
£f) ■ in graphene [8l-tl0|. where it appears in the presence of a random gap in the Dirac spectrum. If the average gap 
04 \ value is small in comparison to the fluctuation strength the system is metallic whereas it is insulating when the gap 
fluctuations are too weak in comparison to the average gap (ill [l2j . 



II. MODEL 



A tight-binding description of electrons in graphene yields the famous energy dispersion with two separate nodes 
r> , (or neutrality points) in the Brillouin zone. In the vicinity of these nodes the momentum dependence of the spectrum 
jj^ ■ is found to be linear and the low-energy behavior of quasi particles can well be described by the Dirac equation 
Hip{x,y) = Eip(x,y) with the Hamiltonian 

H = —ihvp (<r ■ V) + vp m a"3 , (1) 

where vp is the Fermi velocity, a is the vector of Pauli matrices and ip = {<pi, <p-i) is the two component spinor wave 
function, furthermore we set hvp = 1. 

A numerical treatment of the Dirac equation requires a discretization in space. However, the naive discretization 
through replacing the differential operator by a difference operator leads to additional new nodes, which is often called 
fermion doubling or multiplication [Tij| . In real space there are two methods to circumvent this problem 0, Q [H| ■ 
One that we will describe in this section goes back to the idea of Susskind. We start with discretizing the differential 
operator in a symmetric way 

d x f{x) « 2^(//+a - //-a) , (2) 
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where A is the lattice constant which we set to unity in the following. The discrete Dirac equation for m — then 
takes the form 

i i 
~2 ai {^ l+1 > n ~ V'J-i,™} - 2° 2 {^ JLn+1 ~ V'J.n-i} = EaoipLn 

with lattice points given by the coordinates (/, n) with integer I and n. Fourier transformation leads to eigenvalues 
E = ±y/sin(k x ) 2 + siniky) 2 which have four Dirac cones in the Brillouin zone corresponding to four Dirac fermions. 
In order to open a gap at three of them we introduce a lattice operator [l6j which acts on a wave function as 

Blpl,n = ^ {V>M-l,n + tPl-l,n + 1pl,n+l + tpl,n-l} ■ (3) 

Now we discretize Hamiltonian ([T]) by including the lattice operator B and a random gap term 

H^H + 5(B- 2)a 3 + mj,„ a 3 . (4) 

For uniform gap m our new Hamiltonian reads in Fourier representation 

_ fm + 6(cos(k x ) + cos(ky) — 2) sin(k x ) + i sin(k y ) \ 

I sin(fc ;r ) — i sin(fey) —m — 5(cos(k x ) + cos(k y ) — 2)J ^ ' 

with the dispersion 

E = ±^J sin(k x ) 2 + sin(ky) 2 + (m + Scos(k x ) + 8cos(k y ) — 2<5) 2 . (6) 

For m = 0, i5 ^ there is a node at k x = k y — and three additional nodes for m = 0, <5 = at k x , k y = ±ir (cf. Fig.). 
Using this model node degeneracy can be lifted via the parameter S. 

We absorb the index n with the help of matrix representation and write for the wave function 

V'i+i = H Y i[>i + H D ijji-i . (7) 

Each spinor component is now a M-component vector, where M is the width of a strip and thus n = 1, 2, M . The 
matrices H Y , H D read 

Hl„ = 2S" 1 [E a + (26 - m)a 3 ] H Y n+l = S' 1 [ia 2 - 5a 3 ] 

Hn, n -i = -S^ 1 [^2 + 5a 3 ] H® n = -S^ 1 [iai + 6cr 3 ] 

with S = —icf\ + 5a 3 and where H Y has periodic boundary conditions in the y-direction. This matrix structure allows 
us to construct a transfer matrix Ti through the equation [f| 

(^) -*(£)■ « 

The introduction of a different random potential, e.g. random scalar potential, is straight forward. 



A. Lyapunov exponents 

According to 0, 0] the transfer matrices 7], defined in Eq. can be used to calculate Lyapunov characteristic 
exponents (LCE) . With initial values ipo and ipi the iteration of Eq. (|5]) provides i^l by the product matrix 

L 

A/ L =n T '- ( 9 ) 

i=i 

For disordered systems this is a product of random matrices that satisfies Oseledec's theorem (T^]. The latter states 
that there exists a limiting matrix 

r = lim (mIm l ) 1/2L . (10) 

The eigenvalues of T are usually written as exponential functions exp(7i), where 7i is the LCE. Adapting the numerical 
algorithm described in Q, the whole Lyapunov spectrum can be calculated and the smallest LCE is identified with 
the inverse localization length Q. 
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FIG. 1: Brillouin zone of the discrete Dirac equation with circles depicting the positions of the Dirac cones (left). A cut through 
the energy dispersion ((6} is shown on the right for 5 = (blue line), <5 = 0.5 (red line). 



III. NUMERICAL RESULTS FOR RANDOM GAP 



After introducing the model and the corresponding transfer matrices we calculate the inverse of the smallest LCE 
A = l/jmin which is identified as the localization length. A increases with the system width M according to a power 
law: 

A oc M a , (11) 

where a > 1 (a < 1) in the regime of extended (localized) states, and a = 1 in the critical regime. For the exponentially 
localized regime we expect A oc const. According to the one-parameter scaling theory by MacKinnon [l8| . the 
normalized localization length A = A/M obeys the equation 

din A ~. . , 

^lnM=*( lnA )' < 12 > 

where x is an unknown function with solutions of the form 

~A(M,W)=f(Z(W)/M). (13) 

The parameter W characterizes the disorder strength and £ is a characteristic length of the system. The one-parameter 
scaling theory states that A is not depending on M and W separately. Any change of disorder strength W can be 
compensated by a change of the system width M. Furthermore, from the behavior of A in the vicinity of a scale- 
invariant point it is possible to calculate the critical exponent v of the correlation length [f| , which is the localization 
length of the infinite system. This is done by Taylor expansion 



In A = In A r 



Y,A S (\W-W C \M X ' V ) (14) 



s=l 



5 / C \ ~ s l v 



lnA c +E^(i>) ( 15 ) 

s— 1 ^ ' 

with £ = W— W c \~ u . Comparing the latter with eq. (|T3|) . the scaling function £ can be interpreted as the characteristic 
length scale. 



A. Preserved node symmetry: 5 = 



In this case we have a four-fold degeneracy of the node structure. First we calculate A from transfer matrix ([8]) 
with S — 0. If it is not mentioned explicitly we use for the random gap m a box distribution on the interval 
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[rh — W/2, m + W/2], where the corresponding variance is given by W 2 /12. Furthermore, we restrict our calculations 
to the Dirac point (i.e. E = 0). 

Fig. [2] depicts the effect of the average gap m on the localization length A. The localization length always increases 
with system width M, indicating that there is no exponential localization. Only for very weak disorder (W < 0.2) and 
fh = 0.2 the localization length A is almost independent of M, which indicates exponential localization for fh = 0.2 
(cf. Fig. 4(a) |. As disorder increases the localization length decreases monotonically for fh = but not for fh = 0.2 
(cf. Fig. 2(b) If we normalize A by strip width M and perform single parameter scaling as described in [(|, almost 



all data points collapse to a single curve (cf. Fig. 3(a) |. However, we had to neglect data points from weak disorder 
{W < 1.6) to see clearly a scaling behavior. 

The behavior of A for a nonzero average gap (fh = 0.2) is more complex, as shown in Figs. [3(b)| [4(a)] For weak 



disorder the localization length converges to a constant value for increasing M. As disorder increases A increases also 
but remains constant for large M . Then there is a transition at W w 2.1 where A is again growing with system size 
but the slope decreases with increasing disorder. 
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FIG. 2: Localization length for preserved node degeneracy (5 = 0) with average gap m = (left panel) and m = 0.2 (right 
panel) as a function of strip width M. 
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FIG. 3: Scaling plot of the localization length for 5 ■■ 
for W = 0.6; 1.1; 1.6. 



0, m = (left) and 8 = 0, m = 0.2 (right). Left: Rescaled without data 



Due to this behavior of A as a function of disorder it is not possible to perform single parameter scaling in the 
common way. One approach to calculate the scaling function is to minimize the variance of In M — In £ for each 
localization length [18[ . In a double logarithmic plot of A the problem of one parameter scaling translates then into 
shifting all curves onto one @ . Since the position of the resulting curve is irrelevant it is convenient to shift all curves 
onto the lowest i.e. that for biggest disorder. If one looks closely at the data in Fig. 2(b) one sees that this is not 
possible only by shifting. Comparing to Fig. 4(a) one can distinguish two regimes separated at W w 2.1. In both 
regimes one parameter scaling can be performed separately which gives two scaling functions for the infinite system. 
Additionally it is very important to point out that A is always decaying with system size. Usually this is interpreted 
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FIG. 4: Localization length for random gap with zero mean and broken node symmetry with (5 = 0.0, m = 0.2) (left), 5 — 0.5 
and fh — (right) and 8 — 0.5 and m = 0.2 (bottom) as a function of disorder. 

as localizing behavior. Whereas our analysis shows a rather unusual phase transition, namely that the correlation 
length diverges only when approaching the critical point from below W c . In order to extract the functional behavior 
of £ at the transition point wc fitted the data to several functions and found best agreement with 



£(W) oc \W - W c 
The results for the critical parameters are 



for < W < W c ■ 



(16) 



v w 0.289 ±0.013 (W c « 2.156 ±0.009) . 

If we compare the variance g of the fitted critical disorder strength which is g c = 0.387 to the gap width 2m = 0.4 
we see a good agreement. A possible explanation for this might be that if fluctuations of the random gap are larger 
than the gap width states are no more exponentially localized and diffusive transport is possible. From this point of 
view we can also calculate W c from the average gap width which yields W c ~ 2.191. Fitting (|16p with fixed critical 
disorder gives slightly different exponents but also a very good agreement with the numerical scaling function for 
< W < W c : 

v w 0.332 ± 0.004 (W c w 2.191) . 



B. Broken node symmetry: 5^0 



By setting S = 0.5 we break the four-fold degeneracy of the node structure and retain only the node at k x = k y = 0. 
Unlike in the case of 5 = the localization length is not growing with system size if m = 0. Fig. 5(a) shows that for 



weak disorder A is constant with increasing M but decreases with increasing disorder W. However, for W > 4.1 A it 
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increases with M (Fig. |5(b")j ). The normalized data is shown in Fig. |4(b)"| To keep the plot illustrative only a choice 
of the whole data is shown. What can be seen in Fig. |4(b)| is that for weak disorder up to W = 3.6 the normalized 
localization length decays for growing system sizes and scales to zero with M. For disorder larger than W = 3.6 A 
is growing with system size. The growing localization length may be explained by comparison to the clean case. If 
fluctuations of the random gap are in the range of 26 a massless fermion appears. Thus disorder effectively closes the 
gap at the border of the Brillouin zone and the model shows metallic behavior. 
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FIG. 5: Localization length for systems with zero average gap m — and broken node symmetry S = 0.5. Left and right panel 
are for different disorder ranges. 

For weak (i.e. W iS 4) and strong disorder (i.e. W gl 7.5) the behavior is qualitatively the same, characterized by 
a decaying behavior of A with increasing M. The benefit of plotting A over W is that one can see directly two scale 
invariant points where different A are intersecting for all available values of M. These points are indicative of phase 
transitions. Now we use the fitting functions of Eq. (|15)l to extract the critical exponent v from our numerical result. 
For this purpose we set s = 5 and obtain the resulting curves in Fig. [5J The critical parameters are listed in table |U 

Using the scaling form of Eq. (p~3|) all the curves collapse on two curves for a proper choice of the scaling function 
£, as depicted in Figs. [71 [BJ There plots agree qualitatively well for fh = and m = 0.2, the critical exponents for the 
second transition differ slightly though (cf. tables HI and HI)) . 



Critical point 




I 




II 


Exponent v 


1.297 


±0.031 


1.299 


± 0.066 


Wc 


3.975 


± 0.002 


7.668 


± 0.008 


Ac 


0.574 


± 0.009 


0.447 


± 0.005 


Disorder range 


3.87 < 


W < 4.17 


7.35 < 


W < 7.8 


System sizes 


20 < 


M < 80 


30 < 


M < 80 



TABLE I: Critical values for m — and S = 0.5 obtained from fitting the data to equation (|15[) . 



Critical point I II 

Exponent v 1.297 ± 0.045 1.397 ± 0.069 
W c 3.792 ± 0.002 7.629 ± 0.015 

A c 0.591 ±0.007 0.517 ±0.009 

Disorder range 3.72 < W < 3.88 7.1 < W < 8.0 
System sizes 20 < M < 80 20 < M < 80 



TABLE II: Critical values for m = 0.2 and 5 — 0.5 obtained from fitting the data to equation (|15l) . 
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FIG. 7: Rescaled NLL for m = and S — 0.5 around the critical point I (left) and around the critical point II (right). Plots 
contain more data points than used for the fitting procedure to show the validity of one parameter scaling. 



IV. DISCUSSION 



Our numerical results can be summarized as follows. The localization length A always increase with M according 
to the power law of Eq. where the exponent a depends on the model parameters: 



(0 < a < 1 for 5 = 0,m = 

a = for 5 = 0, m = 0.2, W <W C 

< a < 1 for 5 = 0, m = 0.2, W > W c 

a = for S = 0.5, m = 0,0.2, W < W cl 

a>l for 5 = 0.5, m = 0, 0.2, W c x < W < M^ c2 

L0<a<l for 5 = 0.5, m = 0,0.2,11/ > W c2 



(17) 



where <5 = represents the case with four degenerate nodes and 5 = 0.5 a single node. In our numerical results we 
can distinguish these to two cases: (I) For a preserved four-fold node degeneracy (i.e. (5 = 0) the gapless system has 
a monotonically increasing localization length with M as well as with W and does not indicate any transition. In 
the presence of a gap (m ^ 0), however, there is a qualitative change at a characteristic disorder strength W c : For 
W < W c the states are exponentially localized, whereas for W > W c they are not. It is not possible to decide within 
our numerical approach whether they are really extended or power-law localized in the gapped case. As discussed in 
Appendix A, it might be sufficient for diffusion in a 2D system that the states obey a power law. 

(II) For the single node (i.e. S = 0.5) the one-parameter scaling analysis of our results indicates a typical Anderson 
transition at two critical points W c \, W C 2- The exponent a = for weak disorder (i.e. for W < W c \) indicates 
exponentially localized states. There is the intermediate metallic phase for W c \ < W < W C 2 with a > 1 with one- 
parameter scaling behavior near the critical points. This is indicative of two metal-insulator transitions. In particular, 
there is a metal-insulator transition from a = to a > 1 at a critical W c i, which corresponds to a transition from 
Q = 0to0<a<l for the gapped four degenerate Dirac nodes. The difference between a transition from a = 
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FIG. 8: Rescaled NLL for m = 0.2 and S = 0.5 in the vicinity of critical point I (left) and in the vicinity of critical point II 
(right). 

to < a < 1 and a transition from a = to a > 1 is not clear from our numerical results. It could be that the 
latter is a genuine transition from exponentially localized to extended states, whereas the former is a transition from 
exponentially localized to power-law localized states. 

V. CONCLUSION 

We have introduced a model for Dirac fermions on a lattice with several nodes which allows us to perform numerical 
calculations of the localization length within the frame work of the transfer matrix formalism. Using the Hamiltonian 
in Eq. (0 it is possible to break the node symmetry and to compare the properties for one and four nodal points in 
the Brillouin zone. We have shown that states in the gap can be localized and thus the localization length A converges 
to a finite value for increasing system size, whereas in the gapless case there are extended states as expected. 

We have calculated the localization length for various system sizes and for different strength of the random gap. 
In all cases the localization length grows like a power law A ~ M a with increasing system width M. However, the 
exponent a is quite sensitive to the model parameters (cf. (|17[l ). In particular, this exponent vanishes for nonzero 
average gap and weak disorder, indicating exponential localization. Our numerical result also indicates a — for 
non-degenerate nodes, vanishing gap and weak disorder. On the other hand, we have a > 1 only for intermediate 
disorder strength and non-degenerate nodes. Thus, the nodal degeneracy suppresses the intermediate phase. The 
latter is separated from the phases with < a < 1 by transitions that obey one-parameter scaling behavior with 
scale-invariant critical points. This reflects the results of the weak-localization theory, where (anti-)localization has 
been found for (single) two nodes [H, Q . 



Appendix A: Localization and Diffusion in 2D 



Exponentially localized states rule out diffusive behavior. Here we briefly discuss that a power-law decaying state 
can provide diffusive behavior in a 2D electron gas. Diffusion of | 'Q/ (r, | 2 in 2D is defined by the diffusion equation 



d_ D, 

at " 



V 2 |#(r,i)| 2 = 



which has an expanding solution 



|*(r,t)| 2 = K(r,w) 



Dt 



The solution of Eq. (|Aip is also given by the diffusion propagator 

K 



K{q,u) 



u> + Dq 2 



(Al) 



(t - oo) 
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On the other hand, the localization length £ in the spatial direction j can be defined as 



where K(r, uj) is connected with the diffusion propagator by a Fourier transformation: 

Using the Besscl function Jo and the momentum cut-off A for the q integral this result leads to 

K(t,u = 0)-K(t>,u = 0) = §- [ r ^-dx 

D J Xr , x 

and for Ar' , Ar 3> 1 

K (2 f Xr cos(a;-7r/4) 

Thus K(r,uj = 0) decays on large scales like r -1 / 2 . This reflects the fact that a decaying wave function leads to 
diffusion in 2D. 
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